Network disorder and non-afRne deformations in marginal solids 
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The most profound effect of disorder on the elastic response of solids is the non-affinity of local 
displacements whereby the atoms (particles, network junctions) do not simply follow the macroscopic 
strain, as they do in perfect crystals, but undergo additional displacements which result in a softening 
of response. Whether disorder can produce further effects has been an open and difficult question 
due to our poor understanding of non-affinity. Here we present a systematic analysis of this problem 
by allowing both network disorder and lattice coordination to vary continuously under account of 
non-affinity. In one of its limits, our theory, supported by numerical simulations, shows that in 
lattices close to the limit of mechanical stability the elastic response stiffens proportionally to the 
degree of disorder. This result has important implications in a variety of areas: from understanding 
the glass transition problem to the mechanics of biological networks such as the cytoskeleton. 

PACS numbers: 62.20.dc, 83.80.Fg, 81.05Kf 
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I. INTRODUCTION 

In a disordered lattice subject to a small external defor- 
mation, every atom or network junction moves initially in 
an afHne way, that is, displaced in geometric proportion 
to the macroscopic strain. At the same time, the atom 
experiences forces from its nearest-neighbors (NN) which 
are initially displaced afhnely too. In any disordered lat- 
tice atoms will therefore experience a non-zero force due 
to the fact that their NN are randomly oriented. By 
contrast, in a perfectly ordered lattice the forces trans- 
mitted by the NN may balance each other by symmetry 
and the resultant local force in that case can be zero. 
The finite forces acting on each atom in the disordered 
network drive additional local displacements, called non- 
affine displacementsir— . As these are accompanied by 
a release of energy, needed in order to reach mechani- 
cal equilibrium on each atom, the effect of non-affinity 
induced by the disorder on the macroscopic rigidity is 
always to soften the solid with respect to a purely affine 
deformation. This softening is a well established phe- 
nomenon observed many times in simulations and experi- 
mental studies^—. However, since a consistent theoretical 
description of non-affinc clastic response has been elusive 
so far, it has been difficult to predict the rigidity of disor- 
dered solids (including glass, granular media, cytoskeletal 
networks, etc.) as a function of the degree of structural 
order/disorder. Apart from its obvious importance for 
materials science, rigidity is the defining property of the 
solid state of matter— and furthermore the concept of 
(generalized) rigidity as the epitome of cooperative be- 
havior is essential for various other properties, including 
how information is transmitted in many-body systems^. 
In spite of this, the relation between rigidity and disor- 
der has remained largely unknown^. This may well be 
the reason for our present lack of understanding of the 
glass transition phenomenon, wherein rigidity emerges in 
a (disordered) supercooled liquid without any apparent 
change in the symmetry of disordered structure^i^. 

Here we report on a systematic theoretical and numer- 



ical study of the elastic rigidity (described through the 
shear modulus) as a function of the degree of structural 
disorder, which varies from the perfect crystal all the 
way up to the completely disordered lattice. We focus 
our analysis to central-force systems^S where the phys- 
ical bonds between particles have zero bending rigidity, 
and to systems which are close to the vanishing of rigidity 
(marginal .solids). A wide class of physical systems are 
captured by these central force models, such as fiexibly- 
linked biological networks and metallic glasses. In some 
cases, such as covalently bonded materials, junctions can 
possess bending rigidity which makes the transition to 
a rigid phase occur at lower connectivitj*iir— . Without 
loss of generality, our results provide the first direct evi- 
dence for the existence of a broad range of structural dis- 
order where increasing the disorder causes lattices near 
the isostatic limit (i.e. the limit where the number of geo- 
metric constraints exactly matches the number of degrees 
of freedom) to stiffen. This result has deep implications 
in view of the fact that all liquids on their way toward the 
solid state do pass through this isostatic limit, whereby 
even a small difference of mechanical stability with re- 
spect to internal stresses generated by vitrification can 
decide their ultimate fate, i.e. whether to become per- 
manently ordered (crystal) or fully disordered (glass). 



II. THE GAUSSIAN BOND-ORIENTATIONAL 
DISORDER MODEL 



Let us start by introducing the bond-orientation vec- 
tors between nearest-neighbor atoms i and j defined by 
= (cos sin 6', sine/) sin 6*, cos 0). The affine part of 
the shear modulus of harmonic lattices can be derived 
from the expansion of the free energy in the harmonic 
approximation and is given in this limit by the well- 
known Born-Huang formula^, for instance, C^y^y = 
{1/V)kRI J where the lattice sum goes 

over all nearest- neighbor pairs ij. The atomic force 
constant k is defined by the harmonic pair-potential. 



U{rij) = \n{rij — i?o)^, between nearest-neighbors. i?o 
is the interatomic distance at rest in the reference frame 
and V is the volume of the system before deforma- 
tion. Here and below the Roman indices are used to 
label atoms while Greek indices are used to label Carte- 
sian components of vectors, is the connectivity ma- 
trix with Cij = 1 if i and j are nearest neighbors and 
Cij = otherwise in our simple case. Let us consider 
the shear modulus G of a simple cubic lattice in = 3 
in the afhne approximation. Every atom has z = 6 
nearest neighbors and shares a bond with each of them. 
With N atoms in the lattice, on average each atom con- 
tributes z(nf,nj-nf,n^-) to the lattice sum, where the 
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bond-orientation averaged factor is evaluated with an 
appropriate, normalized orientation distribution function 
(ODF) 

= j fi^, 0)sin''6'cos20sin2(/) sin Bdedcj) 

(1) 

Using this averaging the Born-Huang formula can be 
written: 

C^y.y = (2) 

Evaluating f{6, (j)) for different physical situations is one 
of the main exercises of our theory, leading to different 
forms of the resulting shear modulus. For the simple 
cubic system, out of the six orientations of the bonds at- 
tached to each atom, only two are independent, while the 
other 4 can be obtained from these by symmetry opera- 
tions that leave products of terms (a = x, y, z ) un- 
changed. These symmetry operations are either rotations 
(belonging to the octahedral rotation symmetry group 
Oh) about the zenithal (z-)axis, or reflections through 
any plane of symmetry of the cubic cell. One of the two 
independent orientations lies along the zenithal axis and 
is parameterized by O^e = and ^4>ze, while the other 
one lies in the azimuthal plane and is parameterized by 
&az = 7r/2 and 0az = (see Fig. la). The former occurs 
with probability 2/6 and the latter one with probabil- 
ity 4/6. Hence the normalized ODF of the simple cubic 
(SC) lattice can be rewritten, for our purposes, in a sim- 
pler normalized form which is going to be used through- 
out (see the Appendix A for full details): /sc(^, = 
(3/2)[(2/6)J(6l) -I- (4/6)(5(6' - 7r/2)5(</))]. Using this ODF 
in Eq.dl]) immediately leads to {nfjnfjnfjnf^)sc = 
and so the afhne shear modulus for the SC lattice is 
= C^yxy = {NIV)zKRl{n-^nln-^nl)sc = 0. This 
result given by the afhne theory for the perfectly ordered 
simple cubic lattice agrees with Maxwell's counting of 
degrees of frcedoroi^ which indeed predicts G = for 
isostatic lattices with z = 2(i = 6, and also with evi- 
dence from polonium-210 in its alpha form, and the su- 
perconducting calcium-HI, which both feature simple cu- 
bic structure and a vanishing shear stability at ambient 
pressurei^"— . 

In the opposite limit of a completely disordered net- 
work of harmonic springs, one may assume that all bond 
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FIG. 1: (a) Primitive cell of the (perfectly ordered) simple 
cubic lattice. The bonds lying on the horizontal plane are 
called azimuthal bonds (az); those on the vertical axis are 
zenithal (ze). (b) Schematic model of cubic lattice with ran- 
domly added extra bonds (dashed lines), (c) Marginal solids 
with varying degree of disorder parameterized by the degree 
of disorder a defined in the text. 



orientations occur with the same probability. This re- 
sults in a completely de-correlated (bond-orientational) 
disorder described by the uniform ODF fuoiOi 0) = l/47r 
which leads to {nfjnfjnfjnfj)uD = I/I5, where we use the 
subscript UD for "uniform disorder". Hence using this 
resuh in Eq.® gives G = G^^^ = {1/30){N/V)kRIz, 
which is evidently finite for any coordination number 
z. Thus the affine approximation here dramatically fails 
to comply with the vanishing of rigidity at isostaticity 
point (i.e. at z = 6), in contrast with both Maxwell's 
constraint-counting and the empirical evidence^. This 
represents the well-documented failure of the affine ap- 
proximation to describe the elasticity of disordered lat- 

Let us now introduce a way to quantify arbitrary de- 
grees of structural disorder. We again consider the sim- 
ple cubic lattice as our starting point and introduce dis- 
ordered realizations, always keeping z = 6 for the time 
being. The orientations of the bonds now deviate from 
the cubic lattice axes, and we assume that their distribu- 
tion is Gaussian about these axes. A particularly conve- 
nient choice to parameterize Gaussian bond-orientational 
disorder in the network of this topology is the following 
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ODF: fan = (2/6) f^o + 



(4/6) /gl,, with: 



(3) 



where the subscript GD stands for "Gaussian disorder" 
and is the variance, which acts as our quantitative 
measure of lattice disorder. N^eio''^) and Naz{o''^) are 
normahzation factors which also depend on cr^ . We have 
verified that using this compact form for the ODF, in 
terms of the two independent directions only, yields quan- 
titatively the same results that one would obtain by con- 
sidering explicitly all the 6 bond orientations. This ODF 
allows one to span all intermediate degrees of disorder 
by varying the parameter a^. In the limit cr^ — > oo this 
ODF flattens and one recovers the uniform orientational 
distribution /g£)(6', </>) — >• fuD{0,4') = l/47r. In the op- 
posite limit of (T^ — > 0, delta-functions are developed at 
the angles of the cubic lattice so that one recovers the 
simple-cubic ODF, /Gi3(6i,0) fsc{0,(p)- 



III. ANALYTICAL THEORY OF NON-AFFINE 
ELASTIC DEFORMATIONS 

A. Simple cubic and uniformly disordered lattices 

The ultimate cause of non-affinity lies in the fact that 
the particles bonded to a tagged particle, upon the ex- 
ternally applied strain are initially displaced affinely, and 
owing to their displacement they exert a net force on 
the tagged particle. If the particles were placed in an 
ordered fashion around the tagged particle, such as in 
crystals, the resultant sum of these forces would be zero 
due to symmetry. In a disordered solid such a cancelation 
does not occur. Hence the resultant force is finite, and 
it induces an additional displacement on the tagged par- 
ticle in order to keep the mechanical equilibrium, which 
adds to the affine displacement dictated by the imposed 
strain. In this way the mechanical equilibrium is pre- 
served on all particles. Clearly, there is a work associated 
with these non-affine displacements, which bears a neg- 
ative sign since the overall energy is reduced. Formally 
this work represents an additional term in the Born free 
energy expansion and can be written as 



(4) 

The sum runs over all particles and the path integral is 
over the non-affine displacements which the particle 
i undergoes under the action of the force field / exerted 



by the neighbors. H\ 



d^U/drfdr: is the Hessian 



matrix and xf are the coordinates of particle i. Then the 
free energy expansion for a disordered solid accounting 
for non-affinity is: 5F = 5Fa — W where SFa is the 
affine (Born) part of the free energy expansion under a 
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FIG. 2: The theoretical structure-dependent factor in the 
shear modulus, plotted as a function of the degree of Gaussian 
disorder a calculated using Eqs. (1) and (3). The asymptotic 
value of the strong-disorder plateau and the scaling at weak 
disorder are shown by dashed lines. 



given straini^, which leads to the Born-Huang formula 
for the elastic moduli already used above. 

The negative non-affine contribution can be handled 
analytically (further technical details can be found in the 
Appendix A). It has been recently shown by Lemaitre 
and Maloncy by normal mode decompositioi:^ that the 
non-affinc contribution to the modulus takes the form 



r(NA 



k 



V 



(5) 



where S^^ are the affine (force) fields acting on each parti- 
cle, while Vj. and Xk are the eigenvectors (which coincide 
with the normal modes of the solid) and eigenvalues of 
the Hessian matrix, respectively. The sum runs over all 
translational modes (degrees of freedom) of the atoms. 
Both the affine force fields and the eigenvectors y^. de- 
pend on the orientations riy of the bonds in the solid in 
a way geometrically similar to the way the bonds react 
elastically (affinely) to the imposed shear deformation. 
Accordingly, the non-affine term is also proportional to 
the same ("ij'^i'j'^fj"-i'j) factor that has been present in 
the affine model discussed above. As shown in the Ap- 
pendix A, an allowed and convenient form of the eigen- 
modes can be found which leads to the exact evaluation 
of the non-affine terms in the sum in Eq.([5]). Taking the 
average over the disorder of the terms in the sum leads 
to G^^ = C^A- 



(6) 



The important result is that the non-affine part of the 
shear modulus is in fact proportional to the same angu- 
lar average as the affine terms discussed above (because 
the affine fields which are responsible for the non-affine 
term have the same dependence on disorder as the affine 
part), and also to the total number of degrees of freedom 
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in the network, 3A''. which all take part in the force re- 
laxation. The non-afhne effect is, therefore, independent 
of the connectivity z because the non-affine relaxation 
involves the degrees of freedom rather than the bonds. 
In the limit of uniform disorder we may recall that 
= (l/30)(7V/V")Ki?gz. Subtracting the corresponding 
uniform disorder-averaged expression for the non-affine 
part, G^^, we obtain the shear modulus as 



G = KUnZ KUn = 

30 y ° 15 V ° 30V 



kRI{z-Q). (7) 



This equation is derived in the Appendix A. This mod- 
ulus correctly vanishes at z = 6, that is, at the isostatic 
point of central-force solids where the number of con- 
straints (bonds) equals the number of degrees of freedom. 
For z < 6 there is a finite number of floppy modes and 
the systems is no longer rigidi^. 

Equation ([7]) is based on the requirement that the 
Hessian averaged over the orientational disorder be a 
3A^ X 3A^ matrix where the 3x3 submatrices coincide 
with the identity matrix, as shown in the Appendix A. 
This requirement is satisfied for both the simple cubic lat- 
tice and the uniformly disordered lattice that have been 
defined above, i.e. for both the lower and the upper limit 
of the Gaussian disorder spectrum that we introduced in 
the previous section. Even though Eq.© may not hold 
exactly for degrees of disorder intermediate between these 
two limits, we shall continue to assume that Eq.© still 
applies as an approximate relation in between the two 
limits. In the worst case this may give a reasonable in- 
terpolation over the whole disorder spectrum. As this 
is a rather uncontrolled assumption, it will be tested in 
comparison with numerical simulations later in this arti- 
cle. 

Since the same orientation average {n'fjn^jnfjnfj) ap- 
pears as a factor in both G^ and G^^, it follows that 
it plays the same role for an arbitrary degree of disor- 
der, that is, G oc {nfjnfjnfjnfj)cD, where the average 
encodes the information about the lattice disorder via 
fGD{d,(j))- The dependence of the microstructure func- 
tion {nfjnJjnfjn^j)GD upon the degree of disorder as pa- 
rameterized by cr^ is explicitly shown in Fig. 2. In the 
limit CT^ — > oo, the completely uncorrelated uniform dis- 
order case is recovered, {n^j'rv(j'rvljrv(j)GD — >■ 1/15. In the 
opposite limit of weak disorder, the structure function 
goes to zero proportionally to cr^, recovering the simple- 
cubic lattice at (T = 0. The deviation from this propor- 
tionality starts to be significant when ct^ w 1, after which 
the uniform-disorder plateau (= 1/15) sets in. 



B. Gaussian bond-orientational disorder with 
third-body neighbors 

We have therefore established that in the simplest sit- 
uation, where all z NN atoms participate of the same 
degree of disorder, the effects of coordination and disor- 
der are decoupled and they contribute two distinct fac- 



tors to the shear modulus, z — Q (the coordination) and 
(nfjn^^nf^nj'^)G£) (the disorder). An important extension 
arises when bonds between further apart atoms (network 
junctions) are added on a pre-existing lattice. This sit- 
uation describes higher coordination networks and oc- 
curs, for example, in cytoskeletal networks where fila- 
ments are superimposed on an underlying network of a 
different type of filaments. Clearly, the bond orientations 
of the added filaments will be distributed differently from 
the underlying network and thus have a different degree 
of disorder. In order to capture this extension within 
our theory, we consider an underlying isostatic network 
(z = 2d = 6) with a variable degree of disorder onto 
which additional bonds are superimposed and placed at 
random at the network's junctions. This scheme fol- 
lows the most common model employed in the study of 
marginal solida^^^i^ and consists of taking an underlying 
ordered isostatic lattice and then randomly place addi- 
tional springs typically connecting ncxt-nearest-neighbor 
or third-body neighbor atoms (TEN). 

Let us identify the bonds belonging to the underlying 
isostatic lattice as isostatic bonds, while the additional 
bonds will be referred to as excess bonds (see Fig. lb). 
On average, out of z bonds per atom there are 6 iso- 
static bonds of a base cubic lattice, which we assume is 
formed from the variable (Gaussian) disorder, and 2 — 6 
excess bonds per atom. Let us assume the excess bonds 
are placed between TEN atoms along the diagonal of 
the cubic cell, so their average length is of the order of 
\/3i?o- The four third-body neighbor excess bonds in 
the upper half-plane are defined by ^ = 7r/4, V(/), and (f) = 
7r/4, 37r/4, 57r/4, 77r/4 while the four TEN bonds in the 
lower half-plane are parameterized hy 6 = 7r/4-|-7r/2,V(/), 
and (f) ~ 7r/4, 37r/4, 57r/4, 77r/4. Hence, on average a ran- 
domly placed third-body neighbour contributes a factor 
proportional to {nfjTifjnfjn^-) ex ~ 0.06944. This value 
should be affected by the disorder and thus should vary 
with (T. However, it is evident that this value can vary 
only in the very narrow range from 0.06944 (at cr = 0) to 
0.06667 (at a ~ 00). Therefore, the variation of this 
term with a does not lead to a significant change in 
the shear modulus as a function of cr, also because the 
a dependence of the "isostatic" terms is comparatively 
much stronger. This consideration justifies our setting 
{nfjnfjnf^nfj) EX = const = 0.06944 independent of a in 
the following calculations. 



With this decomposition, we can proportionally split 
the affine part of the shear modulus as G^ = G^^ -I- G^ 
where Gf^^ = i{N/V)KRl[6/ z]z{nfjnfjnf^nfj)GD is due 



to the isostatic bonds which occur with frequency 6/z, 
and Gt ^ i{N/V)KRl[{z - 6)/ z]z{nfy^^n-^nl}uD is 
due to the excess bonds that occur with frequency {z — 
6)/z. We can similarly split the non-affine contribution 
into two parts, one due to the 6 isostatic bonds and the 
second due to the additional {z — 6) randomly placed 
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TBN bonds: 



NA 



r-iNA 
^ ex 



= -{N/V)kRIQ 



(8) 



'I] "ij "ij '"tj / GD 



y.) 



^ / X y X y \ 
- \ ^ij n^j n^j n^j ) 



We also need to consider a surface term to account for 
the fact that a finite-size system has the atoms on the 
surface being displaced affinely as they are constrained 
to follow the deformation imposed on the box. The frac- 



tion of atoms on the surface scales with the total number 
of atoms as N~'^/'^ (i.e. sub-extensively). It is diffi- 
cult to precisely determine the number of atoms (the 
thickness of the surface layer) which follow the strain 
of the box in a purely affine way. To account for this 
uncertainty, we introduce a parameter B, of order of 
unity, and write this additional affine surface term as 
G5 = h{N/V)KRlzBN-''/'^{n^jny^n^^jnl)GD- With ah 
these additional contributions put together, the total 
shear modulus takes the form: 



\ ISO 

1 



f^A _ f^NA _ ^ 

ex ISO • 



NA 



(z-6)(n^4n^4)^^ + 3 



(9) 



+ \{N/V)nRlzBN-^/\n-,^nln^^ny^^^ 



GD- 



When the number of atoms N is not infinite, the most 
interesting situation arises for marginal solids as z 6. 
In this limit, and in the limit of a large solid {N ^ 1) 
where we can assume that the surface term is negligible, 
we can also neglect the quadratic term (z — 6)^ and have, 
approximately, G « {N /V)kRI{z - Q){n^jn\-n'^^nlj)GD 
for the non-affine marginal solid with a finite degree of 
disorder measured by a. Guided by the Fig. 2, we can 
split the disorder spectrum into a weak disorder regime 
with < a < 0.1, and a strong disorder regime with 
a > 1. In the strong disorder regime the shear modulus 
is practically constant. Remarkably, in this limit our the- 
ory recovers the result obtained independently from the 
numerical analysis of random packings^! and depleted 
regular latticesi^: G « {l/30){N/V)nRl{z - 6). 

On the other hand, for intermediate and weak disorder 
the shear modulus is strongly dependent upon the degree 
of disorder. For weakly disordered cubic networks, with 
< CT < 0.1, the shear modulus is proportional to a^: 



G « 0.175— Ki?g(z- 6)0-2 (^^Q^ 



If the underlying isostatic lattice is perfectly ordered, i.e. 
CT = 0, the general Ec^. ([9]) has the only contribution aris- 
ing from the excess bonds, giving G = Gex oc (z — 6)^. 
Remarkably, this scaling agrees with the recent result ob- 
tained from the numerical study of square lattices with 
randomly placed extra springs^^. We have now seen that 
the theory predicts that the rigidity of partially disor- 
dered lattices near the isostatic point increases upon in- 
creasing the degree of structural disorder. 



IV. NUMERICAL SIMULATIONS 



To test this conclusion, we performed numerical sim- 
ulations where the disorder is introduced as a perturba- 
tion on a reference cubic lattice (Fig.lc). The magnitude 
of random perturbation for each atom follows the Gaus- 
sian distribution with a variance ct^ along each of the 
three Cartesian axes. The parameter ct, therefore, con- 
trols the degree of quenched disorder in the network, with 
CT = corresponding to the cubic lattice and ct > 1 corre- 
sponding to a completely disordered lattice. The network 
topology is formed by introducing z bonds placed at ran- 
dom on the nearest available atoms iteratively, until each 
atom has z neighbors. A potential energy quadratic in 
the displacement is associated with every bond and the 
total free energy is computed by adding the contributions 
from all the bonds. One should note that the perturba- 
tion of the lattice leads to bond lengths that differ from 
the equilibrium length Rq, an effect which is not consid- 
ered by the theory. However, the contributions from all 
the bonds are summed up and therefore the contributions 
of shorter and longer (than Rq) bonds compensate each 
other yielding something not dissimilar from the equilib- 
rium bond length upon summing up the contributions. 
A further consequence of assuming that all bonds are 
at the equilibrium length is the automatic vanishing of 
the first-order terms in the free-energy Born expansion. 
These terms are related to internal stresses but also in 
this case, upon averaging over the whole network, the 
various contributions cancel to give a small net contri- 
bution (a well-documented fact in the glass simulation 
literature, see e.g. Tanguy et al<^). In weekly connected 
systems (where z < 6 and bond-bending interactions are 
active) the situation might be different as discussed in 
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FIG. 3: The shear modulus of simulated isostatic networks 
{z = 6) as a function of the disorder variance a. Red 
lines: affine, blue lines: non-afRne results. In plot (a) the 
arrows show the curves changing with the increasing net- 
work size, N = 4^,5^,6^,7^,8^,9^ and 10^, respectively. 
In plot (b) the arrows show the convergence of simulation 
(N=216 as an example) for increasing number of steps, t = 
10, 30, 100, 1000, 10^, 310^ and lO''. 



earlier works (e.gji). 

Upon submitting the system to small external defor- 
mations of simple shear, the mechanical equilibrium con- 
dition on each atom/node is enforced such that it fol- 
lows the imposed deformation by moving along paths of 
minimum mechanical energy. This implies, in turn, that 
each atom moves non-afhnely in response to the external 
strain (see Appendix C for detail). From the total free 
energy computed in this way we measure the shear mod- 
ulus of the lattices upon varying the degree of disorder 
a and the connectivity z. The results arc reported in 
Fig. 3. First of all, these plots confirm that our simula- 
tions are valid, both in terms of the scaling G/N ^ Af~2/3 
of Eq.® (Fig. 3a) and the convergence to equilibrium 
(Fig. 3b) tested. It is especially instructive to note that 
when only a short time is allowed for force relaxation, 
the non- affine modulus is very close to the affine result. 

The simulation results reproduce the fundamental laws 
predicted by the theory, as summarized in Eq. ([9]). For 
the isostatic (marginal-solid) limit with z ~ 6, and for 
low to moderate disorder, the shear modulus increases 
in proportion to cr^, and then reaches a strong-disorder 
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FIG. 4: The shear modulus of simulated networks, for a 
fixed value of the disorder variance a — 1, plotted as a 
function of connectivity z, taking only integer values in this 
simulation. The arrow indicates the direction of increasing 
immber of particles: the lattice size from top to bottom is 
TV = 4^ 5^ 6^ 8^ 10^ 12^ 15^ and 20^ The solid line is the 
theoretical scalingi^ G ~ (z — 6), valid for z — > 6. 



plateau at w 1. Further, it is seen that non-affinity 
causes a substantial softening of the solid (lower shear 
modulus) with respect to the purely affine case even 
though the dependence upon a is the same for both 
affine and non-affine simulations. Varying the connec- 
tivity at a fixed degree of disorder, the simulated net- 
works display a rigidity transition at the isostatic point. 
As shown in Fig. 4, the fundamental scaling G cx (z — 6) 
due to non-affinity is asymptotically approached in the 
large limit. In Fig. 5 theory and simulation results are 
shown together as a function of the degree of disorder 
upon varying the connectivity. The theory is in fairly 
good agreement with the simulations for both z = 6 and 
z = 7. It should be noted that the surface term {Gg in 
Eq. ([5])) plays an essential role in this comparison be- 
cause it ensures that the solid remains marginally rigid 
with A^ finite also at z = 6. Further, it is seen that 
upon increasing the connectivity z above the marginal- 
solid limit (z — 6) the shear modulus becomes gradually 
less sensitive to disorder and eventually becomes indepen- 
dent of <7. This is another important and new result: the 
rigidity of central-force lattices with z > 7 is very little 
affected by the degree of structural order/disorder of the 
lattice. In the window of disorder 0.05 < a < 0.5 the the- 
ory appears to underestimate the simulation data which 
exhibit a large hump before reaching the plateau corre- 
sponding to strong disorder. Although there can be many 
reasons for this disagreement, the simulation being com- 
pletely independent, and using a subtly different way to 
define disorder cr, one possibility also lies in the assump- 
tion (Appendix A) that the submatrices of the Hessian 
which determine the form of the eigenvectors (required 
to calculate the non-affine term) reduce to identity ma- 
trices upon averaging, which is an exact procedure only 
for the two opposite SC and UD limits. We have thus 
shown that the rigidity of a cubic lattice increases upon 
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FIG. 5: The shear modulus from the analytical theory (solid 
lines) and simulations (open symbols) upon varying both dis- 
order variance, ct, and lattice connectivity, z. Red symbols 
refer to the affine moduli, blue symbols to the moduli ac- 
counting for non-affinity, for the two discrete cases z = % 
and z = 1 . The continuous variation of z above the isostatic 
point z = 6 shows how the effect of disorder diminishes, and 
gradually disappears at higher connectivity. 



introducing structural bond-oricntational disorder, with 
a law of direct proportionality to the disorder variance. 
Clearly, the disorder-induced stiffening is an effect con- 
trolled by the orientation average (nfjnfjnfjn^j), which 
has a simple geometrical meaning: it measures the extent 
to which the bond-vectors of the lattice are misaligned 
with respect to the shear-deformation and y. In 

fact, the product {nfjnfjnfjnfj) vanishes for the simple 
cubic lattice where all bonds are parallel to either the x or 
the y axis. On the other hand, the product {nfjnfjnfjn^j) 
is maximized and saturates to its highest value (=1/15) 
when the bonds tend to be randomly oriented. In the cu- 
bic lattice the vanishing of (nf.nf.nf.nf,) and the shear 
modulus with it, reflects the fact that there is no restor- 
ing force if all parallel columns in either the x or y di- 
rection are tilted (maintaining the constant separation 
of neighbors), which is an example of a floppy mode. In 
disordered lattices, with bonds misaligned with respect 
to X and y, the longitudinal components of the bond dis- 
placements cannot be neglected (they are no longer sec- 
ond order with respect to the displacement magnitude) 
which implies a finite restoring force and thus a finite 
shear modulus. From a different perspective, the fac- 
tor {nfjTifjnfjn^j) measures the average local (geometric) 
frustration of the bonds and it takes a finite value when- 
ever the bonds no longer can be displaced without paying 
a finite energy cost. 



V. CONCLUSIONS 

First of all, the conclusion that bond-orientational dis- 
order causes stiffening of an ordered lattice can be ex- 
tended to other crystal lattices by re-evaluating the fac- 
tor {nfjn'^jnfjnfj). For example, for the hexagonal close 



packed (hep) lattice with zhcp = 12 wc find that 

'27V3-h3^ 



512 



0.0019. (11) 



As a result, the ratio between the hep modulus and the 
modulus of a uniformly disordered z-coordinated lattice 
takes the form: 



Ghcp _ {'>^'!3n\]n'!jnl^HCP 



G 



0.0285 



Zhcp 



(12) 



where we used that Gna = for the ideal hep lattice. 
If, for instance, the disordered lattice has 2 = 7, then 
Ghcp /GuD ~ 1/3, i.e. the shear modulus of an ideal 
hep lattice with z = 12 (for which the non- affinity is 
certainly negligible) is almost three times smaller than 
the shear modulus of a lattice with uniform disorder and 
z = 7, in spite of the hep lattice being more strongly 
connected. The difference becomes more pronounced 
upon increasing the coordination of the disordered lat- 
tice. Therefore the disorder-induced stiffening of crystal 
lattices is a general effect not limited to the cubic lattice 
studied in detail above. 

Secondly, these findings may have consequences for the 
glass transition of supercooled liquids. Descriptive the- 
ories of the glass transition have emphasized the role of 
local icosahedral symmetry which emerges at the glass 
transition^. Icosahedral clusters have also been identi- 
fied with the cooperatively rearranging regions of mean- 
field (e.g. Adam-Gibbs) theories^i. The growth and jam- 
ming of such icosahedral clusters lead to the sudden in- 
crease of the viscosity and eventually to the structural 
arrest. In the icosahedral clusters an atom at the center 
is bonded to its 12 nearest neighbors, forming an icosa- 
hedron where none of the bonds can be parallel with ei- 
ther the X OT y directions. Instead, clusters with hep 
symmetry could be formed which have the same num- 
ber of bonds as the icosahedra but 6 of these bonds 
can be aligned with the x and y axis. The same ap- 
plies to fee clusters. Therefore it follows that the factor 
{nfjnfjnf^n^j) is much larger for the icosahedra than for 
cither hep or fee clusters and thus the icosahedra are more 
rigid than the hep or fee clusters. This implies that, dur- 
ing the supercooling process, the icosahedra, which are 
the precursor of the rigid glass state^^, are able to stand 
the internal stresses generated by the vitrification^ bet- 
ter than the hep or fee clusters which are the precursor 
of the ordered crystal. This is an unprecedented obser- 
vation that sheds light on the physical origin of the onset 
of rigidity at the glass transition. 

In further developments, this theory can be applied 
to analyze the onset of rigidity in systems which possess 
medium-range order—, such as the technologically impor- 
tant amorphous semiconductors and metallic glasses^SiSi. 
In principle, the orientation distribution function for 
partially disordered lattices (e.g., those containing de- 
fects such as dislocations causing distortion of the bond- 
orientations, or the above mentioned solids with medium- 
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range order) can be extracted from structural investiga- 
tion of the material in question (e.g., scattering or mi- 
croscopy techniques) and used as input in our theory to 
investigate the elastic properties as a function of struc- 
ture. 

In summary, by combining theory and simulations we 
have established a fundamental law which governs the 
rigidity of solids as a function of structural disorder. Ac- 
cording to this law, the stiffness of lattices increases with 
the degree of bond-orientational disorder. Our finding 
suggests a new way of increasing the mechanical sta- 
bility of technologically important marginal-solids such 
as e.g. Po^^° and the superconducting Ca-III, both of 
which exhibit simple cubic lattices which therefore can 
be stabilized by introducing structural disorder (e.g. in 
the form of defects). Furthermore, our findings suggest 
a new way of looking at the unresolved problem of the 
glass transition. In light of our results, vitrifying liquids 
upon crossing the marginal-solidity or isostatic threshold 
find themselves in disordered configurations which are 
intrinsically more rigid than the ordered ones. 
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Appendix A: Non-afRne deformations 

In the following we give a brief outline of the analyti- 
cal method that we devised to calculate the elastic con- 
stants of disordered harmonic lattices which accounts for 
non-affinity2ii^. The method makes use of the non-affine 
linear response formalism for disordered solids which is 
described in detail by Lemaitrc and Malonc}*^. One of 
the central results of the formalism, which represents the 
starting point of our derivation, is the exact lattice sum 
for the elastic moduli accounting for non-affinity^ 

_ RIk , « f X 1 V (^^C • y-k){^^x ■ ^k) 

ij k 

The superscripts A and NA denote the affine and non- 
affine parts of the moduli, respectively. The affine part is 
the standard Born-Huang expression in terms of a lat- 
tice sum over nearest-neighbors (NN) i and j, where 
rijj = (cos (j) sin 9, sin (f) sin 0, cos 0) is the unit vector defin- 
ing the orientation of bond ij . Cij is the adjacency matrix 
(cjj = 1 for two NN atoms and = otherwise) and 
K the harmonic spring constant of the harmonic inter- 
action between two NN atoms. The non-affine part is a 



sum over the 3N eigenmodes of the Hessian or dynamical 
matrix of the lattice. The 3A^ x 3iV Hessian matrix for 
harmonic lattices is given by 

Hf = S., Yl ^^^^<sni - (1 - 6.j)^c.,r^4 (A2) 

s 

Ea. (|A2p follows from replacing the harmonic potential 
U{rij) = \n{rij — Rq)^ in the definition of the Hessian 

matrix: H^j^ = d'^U/drfdr^. Rq is the rest length 
of the bonds, v^. and Afc in Ea. (jAl[) are eigenvectors 
and eigenvalues of the Hessian, respectively. The inner 
product (S^^ • v^) is the projection of the affine force 
field (i.e., the force field exerted on every atom by 
the affine motions of its neighbours) on the eigenvec- 
tor Vj,. The analytical form of the affine fields is given 

by^: S"^^ = ^Z]-^u'*'^u"'S"S"'4j- Thus, the evalua- 

j 

tion of the non-affine term in the elastic moduli reduces 
to the task of evaluating the eigenmodes of the Hessian, 
Xfc ~ Q.r^l^i where r = I...A'^ and I = x,y,z. In gen- 
eral, there are no analytical routes to evaluate the eigen- 
modes. This becomes possible however if one chooses 
mi = 6/ where C; is the standard Cartesian basis of M'^, 
as we show below. This treatment is exact in the case 
of the simple cubic (SC) lattice where the Hessian is ex- 
actly given by H^f = | J2j - (f - %)cy) Sap 

and clearly the Vj. = e; are eigenvectors of the Hes- 
sian because the submatriccs of the Hessian are diagonal. 
In the case of the uniformly disordered (UD) (isotropic) 
lattice as defined in the main article, it can be shown 
that the eigenvectors have the same form. A general the- 
orem states that if v^. and Xk are respectively the eigen- 
vectors and eigenvalues of a matrix A, the same eigen- 
vectors are shared by a matrix f{A), for any function 
/ of the matrix. The eigenvalue equation for the lat- 
ter matrix is then given by: f{A)Y.k = /(^fe)Xfc- Let us 
take /(...) = (...), where (...) denotes the average over 
the bond-orientational disorder, as usual. Then we have 

fiS.) = f [Si-j J2j Cij - (1 - S^j)cij'j 6ai3, where Sap is 

the Kroncckcr's delta. Furthermore, = muj'^ where cok 
are the normal mode frequencies of the solid. The lattice 
with uniform disorder is an isotropic solid, which means 
that the normal mode frequencies depend only on scalar 
quantities. Clearly the same applies to the Afc which in 
turn implies: /(Afc) = Afc. This result establishes that the 
eigenvectors of the form: v^. = C; are also eigenvec- 
tors of the Hessian for the uniformly disordered lattice. 
Hence, in both the SC and UD limits one can write the 
eigenvalue equation for the Hessian as: (ii® /)(a® e) = 

A(a (g) e), where = f (^Sjj J2j ^ij ^ (1 ^ '^y '^"^^^ L 
denotes the 3x3 identity matrix. Then the inner prod- 
ucts of the affine fields with the eigenmodes becomei^: 

arxfc.)= (A3) 

rs,r' s' 
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where the sum runs over two pairs of NN atoms 
at the time, rs and r's'. Consistent with the 
main approximation of this theory, we replace the 
orientation-dependent terms with their orientation av- 
erage. One finds that {nl^n'^gnf^nl,^,n^,^,n^,^,) — 
{Srr'Sss' - 6rs'Ssr')Bi.,^K.x where Bi^,^^^ = {nijU^ijnljnf^) 
and (...) = / ■■■/{&, (f>) s'mddddcf) denotes the average over 
bond orientations according to some orientation distribu- 
tion function f{9,(j)). Replacing this average in Ea.(|A3p. 
one obtains: 



to evaluate the orientation average in Eq. ljASp . In the 
worst situation, this approximation still should provide a 
reasonable interpolation since it correctly describes both 
the lower (SC) and upper (UD) limits of the disorder 
spectrum. The validity of the approximation has been 
checked by comparison with simulation data in the main 
article and as expected it provides a good description of 
the data apart for a window of disorder 0.05 < a < 0.5 
located in the middle of the spectrum. 

Appendix B: Gaussian disorder 



(A4) 



E2 



rs 



Furthermore, we have that 

N 

identities J2'^rJ2'^rs — Yl 



CrsCsr = Crs and thc 

N N 
rs j 



N 



Crs(l — Srs)] = o-rO,sHrs, whcrc wc defined H_ = 

rs 

H^® !_ and / is the 3x3 identity matrix. Recalling 
that J^HrsO-s = Aflr we obtain (S^^ • v^, ) (S^^ • v^, ) = 

s 

KR^XkBi^c^Ky^ Using this result in Eq. ()Aip . the non-affine 
part of thc shear modulus follows as: 

qNA 



V 



k 

N 3 



2 ■ 



(A5) 



V Xk 

k=l 1=1 



which is a key result used in the main part of the article. 

This theory strictly applies to the SC and UD lim- 
its only because only in these specific cases the angular 
average (n^n^) = Sa/sf-i is exact. However, we make 
the additional assumption that for the intermediate de- 
grees of disorder this result still holds approximately, i.e. 
(nfj-nfj.) ~ (5q^/3 for < cr < CX3. Within this approxima- 
tion Ea. (jA5P can be used to evaluate the elastic moduli 
of lattices with disorder (parameterized by variable cr), 
by using the Gaussian ODE defined in the main article 



The distribution function for the orientation of the lat- 
tice bonds (ODE) f{6, (f>) is defined by the following re- 
lation 



..f{9,(t)) sinOdedcj) 



(Bl) 



To evaluate the shear modulus we need to evaluate its 
structure-dependent part 



«j-rif^-n?;,nf^) = j /(6', 0)sin^6'cos>sin^0sin6'd6'd0 

(B2) 

All thc information about the structure of thc lattice is 
contained in J{0,(j}). Eor isostatic lattices in d = 3, as 
explained in the main body of the article, the degree of 
structural disorder can be varied continuously from the 
simple cubic lattice to the complete uncorrelated disorder 
by using the following Gaussian ODE 



az fn 



with 



fa 
JQ 



GD\ 



(B3) 



(B4) 



The normalization factors N^e and Naz are defined by 



/ iV.e(<T')/GD(^, 4>) sinOdedc^ = 1 
J Naz{a^)f --^{9, ^) sin eded^ ^ I 



(B5) 



Combining Eq. (jB3P - (jB5p we obtain the following analyt- 
ical expression for the microstructure-dependent factor 
in thc shear modulus: 
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{nfjn^jnf,n'"^j)GD = J fGD{0,(l))sm'^Ocos'^(l)sm'^(l)s'm0d9d(j) 

= {e-i2'- VT7^v^(Erf[(7r - 5ia'^)/V2^] + Erf[(7r+5icr2)/y2^] 

+5e8'^'(Erf[(7r - 3ia^)/V^] - Erf[(7r + 3ia^)/V^] 

- 2e4'^'(Erf[(7r - ia^)/V2^] - Erf[(7r+icr2)/V2^] + 2iErfi[y^^]) 
+2iErfi[(3Vo^)/V2]) - 2iErfi [5Vg2/V2])} 

X {384(Erf[((7r+icr2)yT72^] - 2iErfi[l/y27^] 
+iErfi[(yT7^(i7r+a2))/V2])}-' 

+ {c-20<T^(5c8'^'(2e4'^'(Erf[(7r - 2ia^)/2V2^] +ETf[{TT+2ia^)/2V2^]) 
+Erf[(^ - 6ia^)/{2V2^)] + ETi[{Tr+6ia^)/{2V2^)]) 
+Erf[(7r - 10ia2)/(2V2^)] 

+Erf[(7r+10ia2)/(2V2^)])( - 2 + 2cS'^'Erf[(V27r)/x/^] 
+Erfc[(y2(^ - 2i(j2))/v^] +Erfc[(y2(7^+2i(T2))/^/^]))} 
X {768Erf[7r/y2^](Erf[(7r - 2icr2)/2^/2^] 
+Erf[(^+2ia2)/(2y2^)])}-' 



Ea. (|B6|) is a function of the Gaussian variance cr^ only 
and has been plotted in Fig. 2 in the main body of the 
article. 



Derivation of the Gaussian ODF 

Here we provide a proof of the expression given by 
Eq.(2) in the main article fEqs.(jB3 |) - (|B4|) in the sec- 
tion above) for the Gaussian ODF used to vary the de- 
gree of disorder in our model. We start by consider- 
ing the orientations of the bonds in the primitive cell of 
the simple cubic (SC) lattice. There are 4 bonds lying 
in the xy plane that are defined by pairs {9, (f>) of an- 
gles (7r/2,0), (7r/2,7r/2), (7r/2,7r), and (7r/2, 37r/2). Fur- 
ther, there are two bonds lying along the polar z-axis 
defined by (O,V0) and (7r,V(/)). The Gaussian model 
realizes the situation where each of these bonds is dis- 
tributed around the simple-cubic angle pairs mentioned 
above according to a Gaussian distribution fan with 
variance a^. In the limit a the model devel- 
ops Dirac deltas at the SC bond orientations thus re- 
covering the perfectly ordered SC lattice. We wish to 
demonstrate that in order to evaluate the angular average 
{nfjn^jn^jn^j)GD = J ...faoid, (j)) smedddcf), one needs 
only to consider the two principal orientations (7r/2,0) 
and (0,V(/)) (from which the other ones can be obtained 
upon application of the symmetry operations of the group 
Oh), which we denote as azimuthal and zenithal respec- 
tively, because all the other orientations contribute terms 
which are identical to either of these two. In other 
words, we want to show that for the purpose of evalu- 
ating (nfjnfjnfjnfj)G£) the following identities hold: 

= g I / . JGD ^ Z^JGD j—QJGD+QJGD 
\P=0 q=0 / 

(B7) 



where we identified: 

/SI, ^ ^ ^,,(a^)e-(''-f )V2.^,-0V2.^ 

fGO^f^T'-^NU^')e-''^'^' 

i.e. the first terms of the two sums in Eq. (jB7p . defined 
by p = and q = 0, respectively. Let us start with the 
orientations in the xy plane. The principal azimuthal 
orientation /S^, with p = 0, contributes to the average 
{nijn\-rifjn\^)GD the following integral: 

j sin'*6lcos20sin20e-(''-^)'/^'"'e-'^'/2<T2 ^^^q^q^^ 

+ J sin'*6lcos^0sin20e-(^~5)'/2'"'e~('^-2'^)'/2<T%jj^5,^5,^0 

^2 J sin'*6lcos20sin2(/)e-(''-5)'/2<T%-0V2a%ij^^^^^^ 

(B9) 

where we omitted the normalization factor and where the 
equality holds because the functions in the two integral 
on the l.h.s. subtend the same area within the same 
interval of integration [0, 27r]. It should be noted that the 
second term on the l.h.s. is strictly required in order to 
have a Gaussian centered on (p ~ 0. Without this term, 
since integration goes from to 2tt, one would have only 
half of the Gaussian and would not properly count the 
contributions from the angles which lie close to (j> — 2tt 
(or equivalently on the negative axis close to = 0~). 

The second term in the first sum in Ea. (|B7[) . with p = 
1, contributes to the average the following term: 

r sm^9e-^'"i^"/^""de r cos2(/,sin20e-(^-^)'/2a'd0 
Jo ^0 

(BIO) 

In this integral, let us change to the variable $ = — -I . 
Upon applying the well-known trigonometric relations: 
sin((/)-t- -1) — cos(t> and cos(0-|- ^) = — sin0, the integral 
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in Eq. (jB10[) becomes: 



(Bll) 



sin^6'e 



This integral is identical to the r.h.s. of Eq. (jB9p because 
the area subtended by the function cos^0sin^(/)e~'^ 1'^" 
(and obviously by cos^$sin^$e^* 1"^" ) in the interval 
[0,27r] is the same as in the interval [— 7r/2, 37r/2], as one 
can easily verify by plotting the functions. This rigor- 
ously establishes the identity /g^ = Jqq^^ = Jq^j ^ ^ 
which holds in the calculation of {n^,rv(,n^rv(^)GD- With 



analogous arguments it is easy to show that Jq}}^^ = 

fcw^""'^^ also holds 
us to conclude that: 



fcD^^^^^ also holds Vp where p is an integer. This leads 



Jgd 

p=0 



4/S 



GD 



(B12) 



To complete the derivation of Ea. (|B7|) (i.e. Eq.(2) 
in the main article), we still need to demonstrate that 

J2q=Q Igd'^'^'^ ~ '^fcD holds as well in the calculation 
of the angular average. The principal zenithal orienta- 
tion foi)-, with g = 0, contributes to the angular average 
{nfjn^jnf^n^j)GD the following integral: 



/ sm^0e~'''/^''^de / cos2(/)sin2 
Jo Jo 



(B13) 



where we omitted again the normalization factor. The 
orientation with q = 1 gives the following contribution 
to the angular average: 



sin'Oe-^'-^^ /"^ d6 I cos^cisin^ 



(B14) 



This integral, upon changing to the variable Q = 9 ~ n, 
becomes: 



/•27r 

-sin^ee-^'/^'^'cie / cos^gisin^ 



where we used sin( 



sind. 



(B15) 



The func- 



tion sin^Oe ® /^"^ is an odd function of the vari- 
able e. This implies that: -sin^ee^^'/^'^'de = 

/J'sin^ee^^'/^'^'de. The last identity estabhshes that 
the integral in Eq. ()B14| ) and the integral in Eq. (|B13p are 
identical. This result in turn demonstrates that, for the 
purpose of calculating {n^n\.rv^.n\.)GD'- 



Er{q-n.y<j>} _ r, fze 
Jgd — ^JGD- 



(B16) 



q=0 



Recollecting the results presented above, we have 
shown rigorously that to the effect of calculating 



{^ii''^ii'"'ii'"'ii)GD the following equality holds: 



Jgd 



igd 



\p=0 



+ l^fGD 

9=0 



£az I rze 

qJgd + qJGD 



6 6 



I ^ AT ^'^2^„-e^/2CT^ 



(B17) 



which has been used in the main article. 

As shown in Fig. 2 of the main article, using Eq.(??) 
to evaluate {nfjnfjnfjnfj)GD leads us to recover exactly 
limc^oc{nfjnJ-nfjnJ-)GD ~ 1/15 which is the value that 
one calculates assuming that the bond orientations are 
randomly distributed in the solid angle, i.e. using the 
uniform-disorder ODE: fjjjj ~ 1/Att. This result gives 
a further confirmation of the correctness of the above 
derivation. 



Appendix C. Simulation of non-affine disordered 
lattices 



We simulate the effect of non-affine deformations on 
the shear modulus of an elastic solid by implementing the 
following scheme. Eirst we place N points on a perfect 
simple cubic lattice, then we perturb the lattice topology 
to induce a variable degree of disorder and we assign a 
harmonic elastic energy to each bond. Einally we deform 
the lattice by implementing a gradient descent method 
at zero-temperature to find the local energy minima. 

Typically we choose a lattice composed of 20 x 20 x 20 
atoms corresponding to A^ = 8000, having confirmed that 
the intrinsically intensive value of the modulus has sat- 
urated on increasing N. Each point is then perturbed 
from its lattice position by drawing three random sam- 
ples from a Gaussian distribution of variance and displac- 
ing the position of each atom by these values in the x, y 
and z directions. This process is repeated independently 
for each of the A^ atoms. We enforce the condition that 
all positions must be in the interval [0, A^] so that all po- 
sitions x,y,z are, in fact, modulo{x, N), modulo{y, N)^ 
and modulo{z, N). Once the positions of the vertices 
have been defined as above, the connections in the net- 
work are introduced by picking a vertex at random and 
forming bonds to each of the z nearest neighbors, pro- 
vided they do not already have z neighbors, i.e. they are 
available. This process is repeated until the network is 
fully formed and no more bonds can be placed. This com- 
pletely defined the topology of the network. The surfaces 
would usually present a problem for this method of form- 
ing connections since at the surface the nearest neighbors 
are predominantly 'behind' the atom in question. We 
overcome this difficulty by employing periodic boundary 
conditions. Eor instance, a vertex in the 'front' y — z sur- 
face, that is, a vertex with a; > can connect to vertices 
in the 'back' y—z surface {x > L) by connecting forwards. 
In essence, instead of there being a single network, there 
is a central network (in which we form connections) , but 
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this network is surrounded by 26 'image networks' which 
connect to the 6 faces, 12 edges and 8 corners such that 
connections from the central network can connect to a 
vertex in any one of the 26 surrounding networks if they 
happen to be one of the z nearest neighbors. In this 
way connections are not biased at the surfaces. Once the 
topology of the network is defined in this way, we deform 
the network affincly by a very small shear strain with 
typically e = lO"'^ (or e = 0.01%). After the afhne de- 
formation we define a non-affine box, which is centered 
on the centre of the original network but extends out to 
90% of the distance in the x,y,z directions compared to 
the original network. Each bond on a vertex exerts a 
force according to Hooke's law. If a vertex lies inside the 
non-affine box, the position is permitted to move under 
the influence of the force in an over-damped way (the 
time-step update is such that the positions are altered 
in the direction of, and with a magnitude proportional 
to, the force). All vertices outside the non-afhne box, 
are not permitted to move and must therefore deform 
affinely. It is via these surface vertices that the strain 
constraint is implemented. This relaxation procedure is 
repeated 104 times and the energy of the network inside 
the non-affine box is recorded as E^^. Each such relax- 
ation procedure is repeated multiple times for each value 
of disorder. The shear modulus of the network is then in- 
ferred from this using the formula: G = 2E^^ /e^ , where 
we have assumed that the contributions to the modulus 
from the terms on the order of can be neglected at 



deformations this small. 

Many of theoretical results in this paper hold in the 
limit of large system size. To verify that the effects ob- 
served in our simulations are not merely finite size effects 
we repeated the simulations over a wide range of system 
sizes {Njnin = 64, Nmax = 8000), and observed that, as 
the system size increases, the curves of modulus against 
disorder approach saturation in a smooth manner, which 
we assume is the limit of infinite system size. Figure 3a 
in the main text illustrates this convergence. 

The method of relaxation used in our simulations nec- 
essarily means that one will never truly reach equilibrium 
in a finite number of time-steps; at best the energy will 
approach its true equilibrium value exponentially. To 
test that the effects observed in this paper were not an 
artefact of incomplete relaxation of local forces, we sim- 
ilarly plotted the curves of modulus against disorder for 
increasing number of relaxation steps (equivalent to the 
time allowed for stress relaxation), ranging from 10 up 
to 10^ in powers of 10^/^. Figure 3b in the main text 
illustrates this convergence. As in the case of finite size 
effects, the effect of finite relaxation time is that upon 
increasing the total number of steps, the curves converge 
smoothly to a limit which we take as the infinite-time 
one. We are therefore confident that neither finite size, 
nor finite time effects are responsible for the observed 
stiffening of network with disorder, but that the effect is 
an inherent property of the network. 
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